border a(t=0,1){x=t; y=-0.5;label=1;}; 
border b(t=0,1){x=1; y=t-0.5;label=1;}; 
border c(t=0,1){x=1-t; y=0.5;label=1;}; 
border d(t=0,1){x=0; y=0.5-t;label=1;}; 

int n=80;
mesh Th=buildmesh (a(n)+b(n)+c(n)+d(n));

fespace Vh(Th,P1);
Vh u,v;

real pi=4*atan(1);
func g = sin(pi*x)*cos(pi*y);
real epsilon = 1e-10;
real mu = 1;
real betax = 100;
real betay = 100;
problem cdr(u,v,solver=UMFPACK) =
    int2d(Th)( mu*u*v+ betax*dx(u)*v+betay*dy(u)*v+epsilon*(dx(u)*dx(v)+dy(u)*dy(v)) )
  - int2d(Th) ( v*(2*pi^2*epsilon*g+betax*pi*cos(pi*x)*cos(pi*y)-betay*pi*sin(pi*x)*sin(pi*y)+ mu*g) ) 
  + on(1,2,3,4,u=g); 

cdr;

cout << "[galerkin] ||u-u_h||_0 = " << sqrt(int2d(Th)( (u-g)^2 )) << "\n";
plot(Th,fill=1,value=1,nbiso=100,wait=1);
plot(u,fill=1,value=1,nbiso=100,wait=1);
v=g;
plot(v,fill=1,value=1,nbiso=100,wait=1);

func delta= 1/(1/hTriangle+epsilon/(hTriangle^2));
//func delta= 1;
problem cdrgals(u,v,solver=UMFPACK) =
    int2d(Th)( mu*u*v+ betax*dx(u)*v+betay*dy(u)*v+epsilon*(dx(u)*dx(v)+dy(u)*dy(v)) )+
    int2d(Th)( delta*(mu*u+betax*dx(u)+betay*dy(u)-epsilon*(dxx(u)+dyy(u)))*
               (mu*v+betax*dx(v)+betay*dy(v)-epsilon*(dxx(v)+dyy(v))) )
  - int2d(Th) ( v*(2*pi^2*epsilon*g+betax*pi*cos(pi*x)*cos(pi*y)-betay*pi*sin(pi*x)*sin(pi*y)+ mu*g) )
  - int2d(Th) ( delta*(mu*v+betax*dx(v)+betay*dy(v)-epsilon*(dxx(v)+dyy(v)))*(2*pi^2*epsilon*g+betax*pi*cos(pi*x)*cos(pi*y)-betay*pi*sin(pi*x)*sin(pi*y)+ mu*g) ) 
  + on(1,2,3,4,u=g); 

cdrgals;

plot(u,fill=1,value=1,nbiso=100,wait=1);
func gx=pi*cos(pi*x)*cos(pi*y);
func gy=pi*sin(pi*x)*sin(pi*y);
func gxx=pi^2*sin(pi*x)*cos(pi*y);
func gyy=pi^2*sin(pi*x)*cos(pi*y);
real normh2 = sqrt(int2d(Th)(g^2+gx^2+gy^2+gxx^2+gyy^2));
cout << "[gals] ||u-u_h||_0 = " << sqrt(int2d(Th)( (u-g)^2 ))/normh2 << "\n";
plot(u,fill=1,value=1,nbiso=100,wait=1);
